1 ReadMe

This is a documentation of the analyses for the manuscript “(Re)Building Trust? Journals’ Open Science Badges Influence Trust in Scientists.” Currenty the data can be obtained from within this file via this link. Once the manuscript is accepted the data will be also available via the OSF.

2 Data import

3 Datawrangling

First we had to reshape data_longest and filter the trustworthiness (METI) and relativism (TSM) items. Then we had to recode the METI items and omit the empty cases.

# A function to inverse the METI items
inv7fct <- function(x) (8-as.numeric(x))

data_psych_prop_METI <- 
data_longest%>%
  filter(substr(variable, 1, 3) == "tru")%>% # see codebook
  mutate(item = substr(variable, 5, 9))%>%
  select(session, item, treat, value)%>%
  pivot_wider(names_from = item,
              values_from = value)%>%
  group_by(session)%>%
  mutate(meas_rep = 1:n(),
         treat = as.factor(treat))%>%
  ungroup()%>%
  mutate_at(vars(exp_1:ben_4),
            list(~inv7fct(.))) # recoding 1-->5, 2-->4, ...))

data_psych_prop_tsm <- 
data_longest%>%
  filter(substr(variable, 1, 3) == "tsm")%>% # see codebook
  mutate(item = substr(variable, 1, 5))%>%
  select(session, item, treat, value)%>%
  pivot_wider(names_from = item,
              values_from = value)%>%
  group_by(session)%>%
  mutate(meas_rep = 1:n(),
         treat = as.factor(treat))%>%
  ungroup()%>%
  mutate_at(vars(starts_with("tsm")),
            list(~as.numeric(as.factor(.))))
  
# Join the data frames
data_psych_prop_joined <- full_join(data_psych_prop_METI, data_psych_prop_tsm)%>%
  mutate(treat = as_factor(treat))
## Joining, by = c("session", "treat", "meas_rep")
# create a list of questionnaire session which are not empty
PID_list <- data_psych_prop_joined%>%
  mutate(frac_NAs = rowMeans(is.na(.)))%>%
  filter(frac_NAs < .8)%>%
  pull(session)%>%
  unique()

data_psych_prop_joined%>%
  mutate(frac_NAs = rowMeans(is.na(.)))%>%
  select(frac_NAs)
## # A tibble: 1,116 x 1
##    frac_NAs
##       <dbl>
##  1    0    
##  2    0    
##  3    0    
##  4    0    
##  5    0.857
##  6    0.857
##  7    0.857
##  8    0.857
##  9    0    
## 10    0    
## # … with 1,106 more rows
# data frame with valid cases (not complete NA)
data_psych_prop <- data_psych_prop_joined%>%
  filter(session %in% PID_list)

# data frames with sum scales
data_scales <- data_psych_prop%>%
  mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
         Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3, 
                                   exp_4, exp_5, exp_6), na.rm = T),
         Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
         Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
         TSM = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), na.rm = T))

data_scales_lables <- data_psych_prop%>%
  mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
         Treatment = fct_recode(Treatment,
                                `Colored badges` = "CB",
                                `Control Condition` = "CC",
                                `Greyed out badges` = "GB"),
         Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3, 
                                   exp_4, exp_5, exp_6), na.rm = T),
         Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
         Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
        `Topic specific multiplism` = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), 
                                               na.rm = T))

data_scales_wide <- data_scales%>%
  select(Experitse, Integrity, Benevolence,
         TSM, Treatment, session)%>%
  gather(Variable, Value, -session, -Treatment)%>%
  mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
  select(-Variable, -Treatment)%>%
  spread(Variable2, Value)

data_scales_lables_wide <- data_scales_lables%>%
  select(Experitse, Integrity, Benevolence,
         `Topic specific multiplism`, Treatment, session)%>%
  gather(Variable, Value, -session, -Treatment)%>%
  mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
  select(-Variable, -Treatment)%>%
  spread(Variable2, Value)

# Data of treatment check
data_tch <- data_longest%>%
  filter(variable %in% c("tch_1", "tch_2", "tch_3", "tch_4"))%>%
  # recoding "weiß nicht[don't know]" as lowest specification of ordinal variable
  mutate(value = as.ordered(value))%>%
  spread(variable, value)%>%
  # remove complete empty cases
  filter(is.na(tch_1) == F & is.na(tch_2) == F &
           is.na(tch_3) == F & is.na(tch_4) == F)%>%
  group_by(session)%>%
  mutate(meas_rep = 1:n())%>%
  ungroup()

data_tch_n <- data_longest%>%
  filter(variable %in% c("tch_1", "tch_2", "tch_3", "tch_4"))%>%
  # recoding "weiß nicht[don't know]" as missing
  mutate(value = as.numeric(ifelse(value == "-999", NA, value)))%>%
  spread(variable, value)%>%
  # remove complete empty cases
  filter(is.na(tch_1) == F & is.na(tch_2) == F &
           is.na(tch_3) == F & is.na(tch_4) == F)%>%
  group_by(session)%>%
  mutate(meas_rep = 1:n())%>%
  ungroup()

4 Sample

4.1 Sample size

length(PID_list)
## [1] 270

4.2 Skipped after first measurement

sum(is.na(data_psych_prop%>%
            filter(meas_rep == 2)%>%
            pull(exp_1)))
## [1] 13

4.3 Demographics

data_longest%>%
  filter(variable %in% c("age", "sex", "semester"))%>%
  spread(variable, value)%>%
  select(-treat)%>%
  mutate(age = as.numeric(age),
         semester = as.numeric(semester),
         sex = as.numeric(sex))%>%
  skim(.)
Data summary
Name Piped data
Number of rows 257
Number of columns 4
_______________________
Column type frequency:
factor 1
numeric 3
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
session 0 1 FALSE 257 TCT: 1, Lw_: 1, NUL: 1, 7z_: 1

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
age 4 0.98 22.89 2.95 18 20 22 25 33 ▇▆▅▁▁
semester 4 0.98 5.86 3.68 1 3 5 9 19 ▇▅▅▁▁
sex 4 0.98 1.71 0.47 1 1 2 2 3 ▃▁▇▁▁

4.3.1 Count on sex

data_longest%>%
  filter(variable %in% c("sex"))%>%
  pull(value)%>%
  table(.)
## .
##   1   2   3 
##  75 176   2

5 Psychometric properties of the measurements

5.1 Descriptive overview

skewn <- function(x) DescTools::Skew(x, na.rm = T)
kurto <- function(x) DescTools::Kurt(x, na.rm = T)
maxabszvalue <- function(x) max(abs(scale(na.omit(x))))

my_skim <- skim_with(numeric = sfl(skewn, kurto, maxabszvalue))

data_scales%>%
  my_skim(.)
Data summary
Name Piped data
Number of rows 540
Number of columns 26
_______________________
Column type frequency:
factor 3
numeric 23
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
session 0 1 FALSE 270 TCT: 2, Lw_: 2, NUL: 2, 7z_: 2
treat 0 1 FALSE 3 GB: 187, CC: 180, CB: 173
Treatment 0 1 FALSE 3 GB: 187, CC: 180, CB: 173

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist skewn kurto maxabszvalue
exp_1 13 0.98 5.43 1.29 1 5.00 6.00 6.0 7 ▁▁▂▃▇ -0.87 0.40 3.42
exp_2 13 0.98 5.44 1.15 2 5.00 6.00 6.0 7 ▂▃▆▇▅ -0.53 -0.21 2.98
exp_3 13 0.98 5.38 1.26 2 5.00 6.00 6.0 7 ▂▃▆▇▅ -0.63 -0.15 2.69
exp_4 13 0.98 5.31 1.32 1 4.00 6.00 6.0 7 ▁▁▂▃▇ -0.74 0.06 3.25
exp_5 13 0.98 5.10 1.39 1 4.00 5.00 6.0 7 ▁▂▃▅▇ -0.60 -0.12 2.96
exp_6 13 0.98 5.43 1.24 1 5.00 6.00 6.0 7 ▁▁▂▃▇ -0.68 0.02 3.57
int_1 13 0.98 5.32 1.24 1 4.00 6.00 6.0 7 ▁▁▃▃▇ -0.49 -0.35 3.48
int_2 13 0.98 5.33 1.31 1 4.00 6.00 6.0 7 ▁▁▃▃▇ -0.47 -0.50 3.30
int_3 13 0.98 5.23 1.17 1 4.00 5.00 6.0 7 ▁▁▅▅▇ -0.10 -0.78 3.60
int_4 13 0.98 5.32 1.22 1 4.00 6.00 6.0 7 ▁▁▃▃▇ -0.45 -0.37 3.53
ben_1 13 0.98 5.23 1.21 1 4.00 5.00 6.0 7 ▁▁▅▅▇ -0.23 -0.45 3.51
ben_2 13 0.98 5.17 1.26 1 4.00 5.00 6.0 7 ▁▁▆▅▇ -0.22 -0.34 3.31
ben_3 13 0.98 5.27 1.26 1 4.00 5.00 6.0 7 ▁▁▃▃▇ -0.63 0.26 3.37
ben_4 13 0.98 5.00 1.23 1 4.00 5.00 6.0 7 ▁▂▆▅▇ -0.09 -0.56 3.26
meas_rep 0 1.00 1.50 0.50 1 1.00 1.50 2.0 2 ▇▁▁▁▇ 0.00 -2.00 1.00
tsm_1 13 0.98 1.87 0.87 1 1.00 2.00 2.0 4 ▇▇▁▃▁ 0.69 -0.34 2.46
tsm_2 13 0.98 2.05 0.95 1 1.00 2.00 3.0 4 ▇▇▁▅▂ 0.51 -0.73 2.06
tsm_3 13 0.98 2.05 0.97 1 1.00 2.00 3.0 4 ▇▇▁▃▂ 0.62 -0.62 2.01
tsm_4 13 0.98 2.29 1.01 1 1.00 2.00 3.0 4 ▆▇▁▆▃ 0.23 -1.06 1.69
Experitse 13 0.98 5.35 1.12 2 4.67 5.50 6.0 7 ▁▂▅▇▅ -0.59 -0.09 3.00
Integrity 13 0.98 5.30 1.08 2 4.50 5.50 6.0 7 ▁▃▆▇▅ -0.28 -0.60 3.05
Benevolence 13 0.98 5.17 1.06 1 4.25 5.25 6.0 7 ▁▁▆▇▆ -0.16 -0.16 3.93
TSM 13 0.98 2.06 0.67 1 1.50 2.00 2.5 4 ▆▆▇▂▁ 0.33 -0.25 2.88

5.2 METI

5.2.1 Dimensionality (CFA)

First we specified two consecutive threedimensional CFA models

# onedimensional model
cfa_meti_model_1d <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                             int_1 + int_2 + int_3 + int_4 +
                             ben_1 + ben_2 + ben_3 + ben_4"

cfa1d_meti_1 <- cfa(cfa_meti_model_1d, data = data_psych_prop%>%filter(meas_rep == 1))
cfa1d_meti_2 <- cfa(cfa_meti_model_1d, data = data_psych_prop%>%filter(meas_rep == 2))

cfa_meti_model_1 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                     int =~ int_1 + int_2 + int_3 + int_4 
                     ben =~ ben_1 + ben_2 + ben_3 + ben_4
                     int_1 ~~ int_2"

cfa_meti_model_2 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                     int =~ int_1 + int_2 + int_3 + int_4 
                     ben =~ ben_1 + ben_2 + ben_3 + ben_4
                     ben_1 ~~ ben_3
                     ben_1 ~~ ben_2"

cfa_meti_1 <- cfa(cfa_meti_model_1, data = data_psych_prop%>%filter(meas_rep == 1))

# define a function which prints the fit
fpf <- function(x){  
  fm_tmp <- fitmeasures(x)
  return(cat(sprintf(
          "χ^2^ = %s, _df_ = %s, CFI = %s, TLI = %s, RMSEA = %s, SRMR = %s, SRMR~between~ = %s, SRMR~within~ = %s",
           round(fm_tmp[c("chisq")],3), 
                 fm_tmp[c("df")],
           round(fm_tmp[c("cfi")],3),
           round(fm_tmp[c("tli")],3),
           round(fm_tmp[c("rmsea")],3),
           round(fm_tmp[c("srmr")],3),
           round(fm_tmp[c("srmr_between")],3),
           round(fm_tmp[c("srmr_within")],3))))
}

# print the fit for cfa_meti_1
fpf(cfa1d_meti_1)

χ2 = 588.932, df = 77, CFI = 0.811, TLI = 0.776, RMSEA = 0.157, SRMR = 0.084, SRMRbetween = NA, SRMRwithin = NA

fpf(cfa_meti_1)

χ2 = 194.174, df = 73, CFI = 0.955, TLI = 0.944, RMSEA = 0.078, SRMR = 0.049, SRMRbetween = NA, SRMRwithin = NA

cfa_meti_2 <- cfa(cfa_meti_model_2, data = data_psych_prop%>%filter(meas_rep == 2))
fpf(cfa1d_meti_2)

χ2 = 936.555, df = 77, CFI = 0.759, TLI = 0.715, RMSEA = 0.208, SRMR = 0.099, SRMRbetween = NA, SRMRwithin = NA

fpf(cfa_meti_2)

χ2 = 212.262, df = 72, CFI = 0.961, TLI = 0.95, RMSEA = 0.087, SRMR = 0.04, SRMRbetween = NA, SRMRwithin = NA

anova(cfa1d_meti_1, cfa_meti_1)
## Chi-Squared Difference Test
## 
##              Df     AIC     BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## cfa_meti_1   73  9738.4  9853.5 194.17                                  
## cfa1d_meti_1 77 10125.1 10225.9 588.93     394.76       4  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cfa1d_meti_2, cfa_meti_2)
## Chi-Squared Difference Test
## 
##              Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## cfa_meti_2   72 8522.8 8639.9 212.26                                  
## cfa1d_meti_2 77 9237.1 9336.5 936.55     724.29       5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In an next step we ran a two-level CFA …

# onedimensional model
mcfa1d_meti_model <- "level: 1
                    meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                            int_1 + int_2 + int_3 + int_4 +
                            ben_1 + ben_2 + ben_3 + ben_4
                    
                    level: 2
                    meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                            int_1 + int_2 + int_3 + int_4 +
                            ben_1 + ben_2 + ben_3 + ben_4"



mcfa_meti_model <- "level: 1
                    exp_w =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                    int_w =~ int_1 + int_2 + int_3 + int_4 
                    ben_w =~ ben_1 + ben_2 + ben_3 + ben_4
                    int_1 ~~ int_2
                    int_3 ~~ ben_4 
                    ben_1 ~~ ben_2 
                 
                    
                    level: 2
                    exp_b =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                    int_b =~ int_1 + int_2 + int_3 + int_4 
                    ben_b =~ ben_1 + ben_2 + ben_3 + ben_4"
mcfa1d_meti <- cfa(mcfa1d_meti_model, data = data_psych_prop, cluster = "session")
mcfa_meti <- cfa(mcfa_meti_model, data = data_psych_prop, cluster = "session")
fpf(mcfa1d_meti)

χ2 = 764.777, df = 154, CFI = 0.894, TLI = 0.875, RMSEA = 0.087, SRMR = 0.271, SRMRbetween = 0.146, SRMRwithin = 0.125

fpf(mcfa_meti)

χ2 = 277.759, df = 145, CFI = 0.977, TLI = 0.971, RMSEA = 0.042, SRMR = 0.172, SRMRbetween = 0.091, SRMRwithin = 0.081

anova(mcfa1d_meti, mcfa_meti)
## Chi-Squared Difference Test
## 
##              Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## mcfa_meti   145 18147 18484 277.76                                  
## mcfa1d_meti 154 18616 18915 764.78     487.02       9  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.2.2 Reliability (McDonalds \(\omega\))

# First Measurement
semTools::reliability(cfa_meti_1)["omega",]
##       exp       int       ben 
## 0.9245381 0.8287494 0.8509624
# Second Measurement
semTools::reliability(cfa_meti_2)["omega",]
##       exp       int       ben 
## 0.9490740 0.9126052 0.8796121

5.3 Topic specific multiplism

5.3.1 Dimensionality (CFA)

First we specified two consecutive onedimensional CFA models

cfa_tsm_model <- "tsm =~ a*tsm_1 + a*tsm_2 + a*tsm_3 + a*tsm_4
                  tsm_2 ~~ tsm_4"

cfa_tsm_1 <- cfa(cfa_tsm_model, data = data_psych_prop%>%filter(meas_rep == 1))
fpf(cfa_tsm_1)

χ2 = 5.302, df = 4, CFI = 0.991, TLI = 0.987, RMSEA = 0.035, SRMR = 0.048, SRMRbetween = NA, SRMRwithin = NA

cfa_tsm_2 <- cfa(cfa_tsm_model, data = data_psych_prop%>%filter(meas_rep == 2))
fpf(cfa_tsm_2)

χ2 = 6.072, df = 4, CFI = 0.99, TLI = 0.985, RMSEA = 0.045, SRMR = 0.055, SRMRbetween = NA, SRMRwithin = NA

In an next step, we ran a two-level CFA …

mcfa_tsm_model <- "level: 1
                    tsm_w =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
                    tsm_2 ~~ tsm_4
                    

                    level: 2
                    tsm_b =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
                    tsm_2 ~~ 0*tsm_2"

mcfa_tsm <- cfa(mcfa_tsm_model, data = data_psych_prop, cluster = "session")
fpf(mcfa_tsm)

χ2 = 4.422, df = 4, CFI = 0.999, TLI = 0.996, RMSEA = 0.014, SRMR = 0.109, SRMRbetween = 0.09, SRMRwithin = 0.019

5.3.2 Reliability (McDonalds \(\omega\))

# First Measurement
semTools::reliability(cfa_tsm_1)["omega",]
## [1] 0.5264177
# Second Measurement
semTools::reliability(cfa_tsm_2)["omega",]
## [1] 0.6459191

5.4 Treatment check

5.4.1 Dimensionality (CFA)

We specified two consecutive onedimensional CFA models

cfa_tch_model <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4
tch_1 ~~ tch_4 "

cfa_tch_1 <- cfa(cfa_tch_model, 
                 data = data_tch_n%>%
                   filter(meas_rep == 1))
fpf(cfa_tch_1)

χ2 = 0.476, df = 1, CFI = 1, TLI = 1.004, RMSEA = 0, SRMR = 0.002, SRMRbetween = NA, SRMRwithin = NA

cfa_tch_2 <- cfa(cfa_tch_model, data = data_tch_n%>%filter(meas_rep == 2))
fpf(cfa_tch_2)

χ2 = 0.055, df = 1, CFI = 1, TLI = 1.018, RMSEA = 0, SRMR = 0.001, SRMRbetween = NA, SRMRwithin = NA

5.4.2 Reliability (Ordinal McDonalds \(\omega\))

# First Measurement
semTools::reliability(cfa_tch_1)["alpha",]
## [1] 0.9444001
# Second Measurement
semTools::reliability(cfa_tch_2)["alpha",]
## [1] 0.9063324

5.5 Table of (M)CFA fit-indices

tibble(`1d CFA METI 1` = fitmeasures(cfa1d_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA METI 2` = fitmeasures(cfa1d_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d CFA METI 1` = fitmeasures(cfa_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d CFA METI 2` = fitmeasures(cfa_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d MCFA METI` = fitmeasures(mcfa1d_meti)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d MCFA METI` = fitmeasures(mcfa_meti)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TSM 1` = fitmeasures(cfa_tsm_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TSM 2` = fitmeasures(cfa_tsm_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d MCFA TSM` = fitmeasures(mcfa_tsm)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TCH 1` = fitmeasures(cfa_tch_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TCH 2` = fitmeasures(cfa_tch_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       rownames = c("χ^2^", "_df_", "CFI", "TLI", "RMSEA", "SRMR", "SRMR~between~", "SRMR~within~", "BIC", "AIC"))%>%
  column_to_rownames(var = "rownames")%>%
  knitr::kable(., digits = 3)
1d CFA METI 1 1d CFA METI 2 3d CFA METI 1 3d CFA METI 2 1d MCFA METI 3d MCFA METI 1d CFA TSM 1 1d CFA TSM 2 1d MCFA TSM 1d CFA TCH 1 1d CFA TCH 2
χ2 588.932 936.555 194.174 212.262 764.777 277.759 5.302 6.072 4.422 0.476 0.055
df 77.000 77.000 73.000 72.000 154.000 145.000 4.000 4.000 4.000 1.000 1.000
CFI 0.811 0.759 0.955 0.961 0.894 0.977 0.991 0.990 0.999 1.000 1.000
TLI 0.776 0.715 0.944 0.950 0.875 0.971 0.987 0.985 0.996 1.004 1.018
RMSEA 0.157 0.208 0.078 0.087 0.087 0.042 0.035 0.045 0.014 0.000 0.000
SRMR 0.084 0.099 0.049 0.040 0.271 0.172 0.048 0.055 0.109 0.002 0.001
SRMRbetween NA NA NA NA 0.146 0.091 NA NA 0.090 NA NA
SRMRwithin NA NA NA NA 0.125 0.081 NA NA 0.019 NA NA
BIC 10225.876 9336.494 9853.512 8639.946 18914.759 18484.147 2791.127 2653.377 5445.587 2079.674 899.974
AIC 10125.120 9237.119 9738.362 8522.826 18616.055 18147.038 2769.537 2632.082 5360.243 2049.381 875.752

6 Results of the treatmentcheck

6.1 Plot

res_tch_data <- data_tch%>%
  gather(variable, value, starts_with("tch_"))%>%
  group_by(treat, variable, value)%>%
  summarize(freq = n())%>%
  ungroup()%>%
  mutate(treat = case_when(treat == "GB" ~ "Greyed out badges",
                           treat == "CB" ~ "Colored badges",
                           treat == "CC" ~ "Control condition",
                           T ~ treat),
         value = case_when(value == "-999" ~ "don't know",
                           T ~ value),
         variable = case_when(variable == "tch_1" ~ "item 1", 
                              variable == "tch_2" ~ "item 2", 
                              variable == "tch_3" ~ "item 3", 
                              variable == "tch_4" ~ "item 4", 
                              variable == "tch_5" ~ "item 5"),
         Frequency = freq)

res_tch_plot <- ggplot(res_tch_data, aes(variable, value)) + 
  geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
  scale_size_continuous(range = c(3,15)) + 
  scale_color_gradient(low = "grey95", high = "grey65") +
  guides(color=guide_legend(), size = guide_legend()) +
  facet_wrap(~treat) + 
  theme_ipsum_ps() + 
  ggtitle("Results of the treatment check", "Frequency per item and experimental condition") + 
  ylab("") + 
  xlab("")
  
res_tch_plot

# res_tch_plot_pub <- ggplot(res_tch_data, aes(variable, value)) + 
#   geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
#   scale_size_continuous(range = c(3,15)) + 
#   scale_color_gradient(low = "grey95", high = "grey65") +
#   guides(color=guide_legend(), size = guide_legend()) +
#   facet_wrap(~treat) + 
#   theme_ipsum_ps() + 
#   ylab("") + 
#   xlab("")

# ggsave("res_tch_plot.svg", res_tch_plot, width = 11, height = 6)
# ggsave("Fig2.jpg", res_tch_plot_pub, width = 110, height = 50, units = "mm", dpi = 300, scale = 2.4)

6.2 Effect size

res_tch_data_A <- data_tch_n%>%
  filter(treat != "CC")%>%
  filter(tch_1 != -999)

effsize::VD.A(tch_1 ~ treat, data = res_tch_data_A)   
## 
## Vargha and Delaney A
## 
## A estimate: 0.8803478 (large)

7 Graphical exploration

7.1 Plot Hyp 1

data_scales_lables%>%
  gather(Variable, Value, Experitse, Integrity, Benevolence)%>%
  ggplot(., aes(x = Treatment, y = Value)) + 
  geom_violin(adjust = 1.5) +
  stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
  coord_flip() +
  facet_wrap(~ Variable, nrow = 1) + 
  labs(title = "Graphical overview (Hyp 1)",
       subtitle = "Violinplots and means ± 1*SD",
       caption = "") +
  ylim(1,7) +
  hrbrthemes::theme_ipsum_ps()
## Warning: Removed 39 rows containing non-finite values (stat_ydensity).
## Warning: Removed 39 rows containing non-finite values (stat_summary).

7.2 Descriptive Effect Sizes Hyp 1

A_GB_CC <- data_scales_lables%>%
  filter(treat != "CB")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate


A_CC_CB <- data_scales_lables%>%
  filter(treat != "GB")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate


A_GB_CB <- data_scales_lables%>%
  filter(treat != "CC")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate
GB CC
CC A = 0.59, d = 0.32
CB A = 0.66, d = 0.57 A = 0.58, d = 0.29

7.3 Hyp 2/3

plot_hyp2_1 <- ggplot(data_scales_lables, 
                      aes(x=`Topic specific multiplism`, y = Integrity)) + 
  geom_jitter() +
  facet_wrap(~ Treatment, nrow = 1) + 
  labs(title = "Graphical overview (Hyp 2/3)",
       subtitle = "Jitter plot per treatment") +
  hrbrthemes::theme_ipsum_ps()


plot_hyp2_1 + stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).

plot_hyp2_1 + stat_smooth(method = "lm")
## Warning: Removed 13 rows containing non-finite values (stat_smooth).

## Warning: Removed 13 rows containing missing values (geom_point).

7.4 Descriptive Effect Sizes Hyp 3/4

Spearman and Kendall correlations:

r_GB <- round(cor(data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_GB <- round(cor(data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

r_CC <- round(cor(data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_CC <- round(cor(data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

r_CB <- round(cor(data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_CB <- round(cor(data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)
GB CC CB
\(r(integrity, multiplism)\) -0.46 -0.36 -0.35
\(\tau(integrity, multiplism)\) -0.33 -0.24 -0.26

7.5 Hyp 4

data_scales_lables%>%
  mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
                               treat == "CB" ~ "Colored badges",
                               treat == "CC" ~ "Control condition",
                               T ~ "treat"))%>%
  ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) + 
  geom_violin(adjust = 1.5, alpha = .5) +
  stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
  coord_flip() +
  labs(title = "Graphical overview (Hyp 4)",
       subtitle = "Violinplots and means ± 1*SD",
       caption = "") +
  xlab("") +
  ylim(1,4) +
  hrbrthemes::theme_ipsum_ps()
## Warning: Removed 13 rows containing non-finite values (stat_ydensity).
## Warning: Removed 13 rows containing non-finite values (stat_summary).

# fig4 <- data_scales_lables%>%
#   mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
#                                treat == "CB" ~ "Colored badges",
#                                treat == "CC" ~ "Control condition",
#                                T ~ "treat"))%>%
#   ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) + 
#   geom_violin(adjust = 1.5, alpha = .5) +
#   stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
#   coord_flip() +
#   xlab("") +
#   ylim(1,4) +
#   hrbrthemes::theme_ipsum_ps()

# ggsave("Fig4.jpg", fig4, width = 120, height = 70, units = "mm", dpi = 300, scale = 1.5)

7.6 Descriptive Effect Sizes Hyp 4

A_mult_GB_CC <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "CB")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_GB_CC <- qnorm(A_mult_GB_CC)*sqrt(2)


A_mult_CC_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "GB")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_CC_CB <- qnorm(A_mult_CC_CB)*sqrt(2)


A_mult_GB_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "CC")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_GB_CB <- qnorm(A_mult_GB_CB)*sqrt(2)
GB CC
CC A = 0.43, d = -0.26
CB A = 0.43, d = -0.25 A = 0.5, d = 0.01

8 Inference statistics (Hyp 1)

8.1 Description of the variables

All of the following analyses are based on the data frame object data_scales_wide which is why we describe it here somewhat more detailed.
All analyses are based on measured variables Integrity (Integrity, is information source sincere, honest, just, unselfish and fair?) and Topic Specific Multiplism (TSM, is knowledge and knowing about a topic arbitrary?). As this data set is in wide format, the experimental conditions are encoded wtihin the variable names:

  • GB means, that the participants of the study could learn from the grey out badges (ans corresponding explanations) within the abstract, that the authors of the study denied to use Open Practices
  • CC means, that the participants of the study could not learn if or if not the authors of the study used Open Practices
  • CBmeans, that the participants of the study could learn from the colored badges (ans corresponding explanations) within the abstract, that the authors of the study used Open Practices

Finally, the variable session identified the study participants.

If we look descriptively at these variables:

data_scales_wide%>%
  my_skim(.)
Data summary
Name Piped data
Number of rows 270
Number of columns 13
_______________________
Column type frequency:
factor 1
numeric 12
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
session 0 1 FALSE 270 TCT: 1, Lw_: 1, NUL: 1, 7z_: 1

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist skewn kurto maxabszvalue
Benevolence_CB 100 0.63 5.44 0.99 2.75 4.75 5.25 6.25 7 ▁▂▇▅▆ -0.08 -0.71 2.71
Benevolence_CC 97 0.64 5.20 0.96 2.75 4.50 5.25 6.00 7 ▁▅▇▆▃ 0.03 -0.61 2.55
Benevolence_GB 86 0.68 4.89 1.15 1.00 4.00 4.75 5.75 7 ▁▁▇▆▅ -0.14 -0.02 3.39
Experitse_CB 100 0.63 5.63 1.10 2.33 4.88 5.83 6.50 7 ▁▂▃▇▇ -0.78 0.13 2.99
Experitse_CC 97 0.64 5.35 0.97 2.50 4.67 5.33 6.00 7 ▁▃▆▇▅ -0.25 -0.48 2.94
Experitse_GB 86 0.68 5.09 1.19 2.00 4.29 5.33 6.00 7 ▂▂▅▇▃ -0.59 -0.32 2.59
Integrity_CB 100 0.63 5.62 0.98 3.25 5.00 5.75 6.50 7 ▂▅▇▇▇ -0.31 -0.83 2.42
Integrity_CC 97 0.64 5.34 0.97 2.75 4.50 5.50 6.00 7 ▁▃▇▇▅ -0.22 -0.62 2.66
Integrity_GB 86 0.68 4.97 1.18 2.00 4.00 5.00 6.00 7 ▂▅▆▇▃ -0.09 -0.76 2.53
TSM_CB 100 0.63 2.01 0.67 1.00 1.50 2.00 2.50 4 ▇▇▇▂▁ 0.38 -0.34 2.95
TSM_CC 97 0.64 1.99 0.61 1.00 1.50 2.00 2.50 4 ▆▆▇▁▁ 0.25 -0.12 3.30
TSM_GB 86 0.68 2.18 0.72 1.00 1.75 2.25 2.56 4 ▅▆▇▂▂ 0.24 -0.46 2.54

8.2 Investigating the missingness

8.2.1 Missingness per Variable

library(naniar)
## 
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
## 
##     n_complete
visdat::vis_miss(data_scales_wide%>%select(-session)) + 
  ggtitle("Missingness per Variable") + 
  theme(plot.margin = margin(1, 2, 1, 1, "cm"))

8.2.2 Marginal distributions Integrity_GB

library(patchwork)
ggplot(data_scales_wide, aes(Integrity_GB, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

8.2.3 Marginal distributions Integrity_CC

ggplot(data_scales_wide, aes(Integrity_GB, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

8.2.4 Marginal distributions Integrity_CB

ggplot(data_scales_wide, aes(Integrity_GB, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

8.2.5 Marginal distributions TSM_GB

ggplot(data_scales_wide, aes(Integrity_GB, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_GB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

8.2.6 Marginal distributions TSM_CC

ggplot(data_scales_wide, aes(Integrity_GB, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_CC)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

8.2.7 Marginal distributions TSM_CB

ggplot(data_scales_wide, aes(Integrity_GB, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_CB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

8.2.8 Cohen’s d of missing/not-missing per variable

d_matrix_missings <- miceadds::mi_dstat(data_scales_wide%>%select(starts_with("Int"), starts_with("TSM")))%>%
  round(.,4)

knitr::kable(d_matrix_missings, caption = "Boxplot of Cohen's d of missing/not-missing per variable")
Boxplot of Cohen’s d of missing/not-missing per variable
Integrity_CB Integrity_CC Integrity_GB TSM_CB TSM_CC TSM_GB
Integrity_CB NaN -0.2103 -0.0096 NaN -0.1009 -0.3312
Integrity_CC -0.1596 NaN -0.0074 0.1056 NaN 0.3622
Integrity_GB 0.1045 0.1696 NaN -0.1339 0.1018 NaN
TSM_CB NaN -0.2103 -0.0096 NaN -0.1009 -0.3312
TSM_CC -0.1596 NaN -0.0074 0.1056 NaN 0.3622
TSM_GB 0.1045 0.1696 NaN -0.1339 0.1018 NaN 
boxplot(as.vector(d_matrix_missings))
title("Boxplot of Cohen's d of missing/not-missing per variable")

8.3 Imputation

M <- 1000
out <- mice(data = data_scales_wide%>%select(-session),
            m = M,
            meth=c("norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm"), 
            diagnostics = FALSE,
            printFlag = F,
            seed = 83851)

8.4 Check of first 10 imputation

out_first10 <-   mice(data = data_scales_wide%>%select(-session),
            m = 10,
            meth=c("norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm"), 
            diagnostics = FALSE,
            printFlag = F,
            seed = 83851)

densityplot(out_first10)

plot(out_first10)

8.5 Parameter and BF estimation

# Set up the matrices for the estimates ##############
# setup of matrices to store multiple estimates
mulest_hyp1 <- matrix(0,nrow=M,ncol=3) 
# and covariance matrices
covwithin_hyp1 <- matrix(0,nrow=3,ncol=3) 


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 within_hyp1 <- lm(cbind(Integrity_GB,Integrity_CC,Integrity_CB)~1, 
                   data=mice::complete(out,i)) # estimate the means of the three variables
 mulest_hyp1[i,]<-coef(within_hyp1)[1:3] # store these means in the matrix `mulres`
 covwithin_hyp1<-covwithin_hyp1 + 1/M * vcov(within_hyp1)[1:3,1:3] # compute the averages 
}


# Compute the average of the estimates ###############
estimates_hyp1 <- colMeans(mulest_hyp1)
names(estimates_hyp1) <- c("Integrity_GB","Integrity_CC","Integrity_CB")
covbetween_hyp1 <- cov(mulest_hyp1) # between covariance matrix
covariance_hyp1 <- covwithin_hyp1 + (1+1/M)*covbetween_hyp1 # total variance


# Determine the effective and real sample sizes ######
samp_hyp1 <- nrow(data_scales_wide) # real sample size
nucom_hyp1<-samp_hyp1-length(estimates_hyp1)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp1 <- (1+1/M)*(1/length(estimates_hyp1))* 
  sum(diag(covbetween_hyp1 %*% 
             MASS::ginv(covariance_hyp1))) # ... (43)
nuold_hyp1<-(M-1)/(lam_hyp1^2) # ... (44)
nuobs_hyp1<-(nucom_hyp1+1)/(nucom_hyp1+3)*nucom_hyp1*(1-lam_hyp1) # ... (46)
nu_hyp1<- nuold_hyp1*nuobs_hyp1/(nuold_hyp1+nuobs_hyp1) # ... (47)
fracmis_hyp1 <- (nu_hyp1+1)/(nu_hyp1+3)*lam_hyp1 + 2/(nu_hyp1+3) # ... (48)
neff_hyp1<-samp_hyp1-samp_hyp1*fracmis_hyp1 # = 172 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp1<-list(covariance_hyp1)


# Test the hypotheses with bain ######################
results_hyp1 <- bain(estimates_hyp1,
               "Integrity_GB<Integrity_CC<Integrity_CB;
                Integrity_GB=Integrity_CC=Integrity_CB;
                Integrity_GB<Integrity_CC=Integrity_CB",
                n = neff_hyp1, Sigma=covariance_hyp1,    
                group_parameters=3, joint_parameters = 0)
print(results_hyp1)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit   Com   BF.u  BF.c     PMPa  PMPb 
## H1 0.999 0.178 5.608 4426.129 0.974 0.830
## H2 0.000 0.279 0.000 0.000    0.000 0.000
## H3 0.039 0.259 0.152 0.152    0.026 0.022
## Hu                                  0.148
## 
## Hypotheses:
##   H1: Integrity_GB<Integrity_CC<Integrity_CB
##   H2: Integrity_GB=Integrity_CC=Integrity_CB
##   H3: Integrity_GB<Integrity_CC=Integrity_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
summary(results_hyp1)
##      Parameter        n Estimate       lb       ub
## 1 Integrity_GB 169.8074 4.992651 4.829868 5.155435
## 2 Integrity_CC 169.8074 5.326883 5.189787 5.463978
## 3 Integrity_CB 169.8074 5.586284 5.444758 5.727809
results_hyp1$BFmatrix
##              H1       H2           H3
## H1 1.000000e+00 49594808 3.700934e+01
## H2 2.016340e-08        1 7.462341e-07
## H3 2.702021e-02  1340062 1.000000e+00

9 Inference statistics (Hyp 2)

9.1 Parameter estimation

path_mod <- "Integrity_GB ~ TSM_GB
             Integrity_CC ~ TSM_CC             
             Integrity_CB ~ TSM_CB"

# Set up the matrices for the estimates ##############
best_hyp2 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp2 <- matrix(0,nrow=3,ncol=3) # and covariance matrices


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 path_fit <- sem(path_mod, 
                 data = mice::complete(out, i), 
                 std.lv = T) # estimate the path coefficients 
 best_hyp2[i,] <- parameterestimates(path_fit, standardized = T)%>% # store path coefficients
                     filter(op == "~")%>%
                     pull(std.all) 
 covwithin_hyp2 <- covwithin_hyp2 + # compute the average of the covariance matrices
                  1/M * lavInspect(path_fit, 
                                   "vcov.std.all")[c("Integrity_GB~TSM_GB", 
                                                     "Integrity_CC~TSM_CC", 
                                                     "Integrity_CB~TSM_CB"),
                                                   c("Integrity_GB~TSM_GB", 
                                                     "Integrity_CC~TSM_CC", 
                                                     "Integrity_CB~TSM_CB")]
}

# Compute the average of the estimates ###############
estimates_hyp2 <- colMeans(best_hyp2)
names(estimates_hyp2) <- c("Int_on_TSM_GB", "Int_on_TSM_CC", "Int_on_TSM_CB")
round(estimates_hyp2, 2)
## Int_on_TSM_GB Int_on_TSM_CC Int_on_TSM_CB 
##         -0.44         -0.41         -0.32

9.2 Visual path model

covbetween_hyp2 <- cov(best_hyp2) # between covariance matrix
covariance_hyp2 <- covwithin_hyp2 + (1+1/M)*covbetween_hyp2 # total variance

# Determine the effective and real sample sizes ######
samp_hyp2 <- nrow(data_scales_wide) # real sample size
nucom_hyp2 <- samp_hyp2-length(estimates_hyp2)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp2 <- (1+1/M)*(1/length(estimates_hyp2))* 
  sum(diag(covbetween_hyp2 %*% 
             MASS::ginv(covariance_hyp2))) # ... (43)
nuold_hyp2 <- (M-1)/(lam_hyp2^2) # ... (44)
nuobs_hyp2 <- (nucom_hyp2+1)/(nucom_hyp2+3)*nucom_hyp2*(1-lam_hyp2) # ... (46)
nu_hyp2 <- nuold_hyp2*nuobs_hyp2/(nuold_hyp2+nuobs_hyp2) # ... (47)
fracmis_hyp2 <- (nu_hyp2+1)/(nu_hyp2+3)*lam_hyp2 + 2/(nu_hyp2+3) # ... (48)
neff_hyp2 <- samp_hyp2-samp_hyp2*fracmis_hyp2 # = 114 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp2 <- list(covariance_hyp2)

9.3 Bayes Factor estimation (Hyp 2)

results_hyp2 <- bain(estimates_hyp2, 
                     "Int_on_TSM_GB < 0 & Int_on_TSM_CC < 0 & Int_on_TSM_CB < 0;
                     Int_on_TSM_GB = 0 & Int_on_TSM_CC = 0 & Int_on_TSM_CB = 0",
                    Sigma = covariance_hyp2,
                    n = neff_hyp2,
                    group_parameters = 3,
                    joint_parameters = 0,
                    standardize = F)
print(results_hyp2)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit   Com   BF.u  BF.c         PMPa  PMPb 
## H1 1.000 0.158 6.326 12953840.149 1.000 0.864
## H2 0.000 1.285 0.000 0.000        0.000 0.000
## Hu                                      0.136
## 
## Hypotheses:
##   H1: Int_on_TSM_GB<0&Int_on_TSM_CC<0&Int_on_TSM_CB<0
##   H2: Int_on_TSM_GB=0&Int_on_TSM_CC=0&Int_on_TSM_CB=0
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp2$BFmatrix
##             H1           H2
## H1 1.00000e+00 8.214022e+21
## H2 1.21743e-22 1.000000e+00

10 Inference statistics (Hyp 3)

results_hyp3 <- bain(estimates_hyp2, 
                    "Int_on_TSM_GB < Int_on_TSM_CC < Int_on_TSM_CB;
                     (Int_on_TSM_GB, Int_on_TSM_CC) < Int_on_TSM_CB;
                     Int_on_TSM_GB = Int_on_TSM_CC = Int_on_TSM_CB",
                    Sigma = covariance_hyp2,
                    n = neff_hyp2,
                    group_parameters = 3,
                    joint_parameters = 0,
                    standardize = T)

print(results_hyp3)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit    Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.541  0.177 3.058  5.488  0.122 0.118
## H2 0.820  0.354 2.319  8.332  0.093 0.089
## H3 10.086 0.513 19.645 19.645 0.785 0.755
## Hu                                  0.038
## 
## Hypotheses:
##   H1: Int_on_TSM_GB<Int_on_TSM_CC<Int_on_TSM_CB
##   H2: (Int_on_TSM_GB,Int_on_TSM_CC)<Int_on_TSM_CB
##   H3: Int_on_TSM_GB=Int_on_TSM_CC=Int_on_TSM_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp3$BFmatrix
##           H1       H2        H3
## H1 1.0000000 1.318921 0.1556656
## H2 0.7581957 1.000000 0.1180250
## H3 6.4240278 8.472784 1.0000000

11 Inference statistics (Hyp 4)

# Set up the matrices for the estimates ##############
mulest_hyp4 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp4 <- matrix(0,nrow=3,ncol=3) # and covariance matrices


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 within_hyp4 <- lm(cbind( TSM_GB, TSM_CC, TSM_CB) ~ 1, 
                   data=mice::complete(out,i)) # estimate the means of the three variables
 mulest_hyp4[i,]<-coef(within_hyp4)[1:3] # store these means in the matrix `mulres`
 covwithin_hyp4<-covwithin_hyp4 + 1/M * vcov(within_hyp4)[1:3,1:3] # compute the averages 
}


# Compute the average of the estimates ###############
estimates_hyp4 <- colMeans(mulest_hyp4)
names(estimates_hyp4) <- c("TSM_GB","TSM_CC","TSM_CB")
covbetween_hyp4 <- cov(mulest_hyp4) # between covariance matrix
covariance_hyp4 <- covwithin_hyp4 + (1+1/M)*covbetween_hyp4 # total variance


# Determine the effective and real sample sizes ######
samp_hyp4 <- nrow(data_scales_wide) # real sample size
nucom_hyp4<-samp_hyp4-length(estimates_hyp4)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp4 <- (1+1/M)*(1/length(estimates_hyp4))* 
  sum(diag(covbetween_hyp4 %*% 
             MASS::ginv(covariance_hyp4))) # ... (43)
nuold_hyp4<-(M-1)/(lam_hyp4^2) # ... (44)
nuobs_hyp4<-(nucom_hyp4+1)/(nucom_hyp4+3)*nucom_hyp4*(1-lam_hyp4) # ... (46)
nu_hyp4<- nuold_hyp4*nuobs_hyp4/(nuold_hyp4+nuobs_hyp4) # ... (47)
fracmis_hyp4 <- (nu_hyp4+1)/(nu_hyp4+3)*lam_hyp4 + 2/(nu_hyp4+3) # ... (48)
neff_hyp4<-samp_hyp4-samp_hyp4*fracmis_hyp4 # = 172 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp4<-list(covariance_hyp4)


# Test the hypotheses with bain ######################
results_hyp4 <- bain(estimates_hyp4,
               "TSM_GB=TSM_CC=TSM_CB;
                TSM_GB>TSM_CC>TSM_CB;
                (TSM_GB,TSM_CC)>TSM_CB",
                n = neff_hyp4, Sigma=covariance_hyp4,    
                group_parameters=3, joint_parameters = 0)

print(results_hyp4)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit   Com   BF.u  BF.c  PMPa  PMPb 
## H1 0.293 0.507 0.578 0.578 0.094 0.081
## H2 0.657 0.182 3.609 8.598 0.589 0.507
## H3 0.657 0.339 1.938 3.734 0.316 0.272
## Hu                               0.140
## 
## Hypotheses:
##   H1: TSM_GB=TSM_CC=TSM_CB
##   H2: TSM_GB>TSM_CC>TSM_CB
##   H3: (TSM_GB,TSM_CC)>TSM_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp4$BFmatrix
##          H1        H2        H3
## H1 1.000000 0.1602312 0.2984187
## H2 6.240981 1.0000000 1.8624253
## H3 3.350996 0.5369343 1.0000000